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Calculations of the photonic band structure, transmission coefficients, and quality factors of various two- 
dimensional, periodic and aperiodic, dielectric photonic crystals by using the finite element method (FEM) are 
reported. The fundamental equations governing the propagation of electromagnetic waves in inhomogeneous 
media are revisited together with the boundary conditions required for each of the performed calculations. A 
detailed account of the eigenvalue and harmonic propagation analysis of the electromagnetic problem is reported 
for several periodic and finite-length structures. It is found that this method reproduces quite well previous 
results for these lattices obtained with the standard plane wave method with regards to the eigenvalue analysis 
(photonic band structure calculations). However, in contrast with frequency methods, the finite element method 
easily allows one to study the time-harmonic propagation of electromagnetic fields and, thus, to calculate the 
transmission coefficient of finite clusters in a natural way. Moreover, the advantages of using this real space 
method for structures of arbitrary complexity are also discussed. In addition, point defect cluster quality factor 
calculations are reported by means of FEM and they are compared to the ones obtained with the FDTD and 
Harminv methods. As a result, FEM comes out as an effective, stable, robust, and rigorous tool to study light 
propagation and confinement in both periodic and aperiodic dielectric photonic crystals. 



I. INTRODUCTION 

Photonic crystals have generated a surge of interest in the 
last decades because they offer the possibility to control the 
propagation of light to an unprecedented level' ^. In its sim- 
plest form, a photonic crystal is an engineered inhomogeneous 
periodic structure made of two or more materials with very 
different dielectric constants. When an electromagnetic wave 
(EM) propagates in such a structure whose period is compa- 
rable to the wavelength of the wave, unexpected behaviors 
occur Among the most interesting ones are the possibility 
of forming a complete photonic band gap (CPBG)^ ^, that is, 
a frequency range for which no photons having frequencies 
within that range can propagate through the photonic crystal 
(PC), to localize light by introducing several types of defects 
in the lattice, or enhancing certain non-linear phenomena due 
to small or anomalous group velocity effects. Unfortunately, 
all these nice features come at a price: the length scales re- 
quired in order to fabricate a photonic crystal appropriate for 
operation at telecommunication frequencies are below the mi- 
cron, so that ingenious innovations were required in order to 
actually fabricate such structures. Furthermore, fabricating 
devices based on these lattices is even more challenging from 
a technological point of view, especially in three dimensions 
. This has resulted in a lot of effort being devoted to investi- 
gating photonic crystal devices based on two-dimensional het- 
erostructures, such us three-dimensional waveguides made of 
a two-dimensional photonic crystal core sandwiched between 
two layers of substrate that confines light by simple refraction 
index matching. Therefore, investigating two-dimensional 
photonic crystals is not a mere academic exercise but an im- 
portant task, both at the fundamental and applied levels. From 
a theoretical perspective, one of the most important venues of 
research is to develop numerical methods to solve Maxwell's 
equations of an EM wave propagating in a PC that are reli- 
able, fast, and capable of dealing with large systems as close 



as possible to the real ones employed for experiments. The 
reason why numerical methods are so important in this field is 
rooted on the fact that the predictions of Maxwell's equations 
for these systems are in excellent qualitative and quantitative 
agreement with experiments and they are able to describe the 
ample phenomenology exhibited by these systems. The de- 
velopment of these numerical methods has undergone an evo- 
lution that mimics to a certain extent the one underwent by 
the techniques used to calculate the electronic band structure 
of semiconductors. Thereby, while early research efforts fo- 
cused on ideal infinite periodic photonic crystals, for which 
techniques on the frequency space are ideally suited (such as 
the plane wave method (PW)23i or algorithms based on the 
tight binding method), it was soon realized that the real in- 
terest of these systems is disordered photonic crystals. Fre- 
quency methods can be used to address the calculation of cer- 
tain quantities (such as the photonic band structure or density 
of states) for some types of disorder, mainly point or line de- 
fects, by performing supercell calculations'* but they are not 
very useful for studying disorder that is not localized (such 
as random displacements of the dielectric constituents of the 
photonic crystal or compositional disorder), which are of great 
practical interest, as these types of disorder are usually asso- 
ciated to the fabrication techniques themselves. Furthermore, 
frequency methods, in principle, are only applicable to infinite 
systems whereas real systems are always finite. For these rea- 
sons, people have started to pay attention during the last years 
to real space methods. These methods diff'er from the previ- 
ous ones in that they work with a finite cluster which does not 
need to be periodic. 

Interestingly, a method that has received very little atten- 
tion in this field even though is has been known and exten- 
sively used in other areas of Physics and Engineering is the 
finite element (EE) method'*' Indeed, there are a number 
of reasons that seem to suggest that it could be applied well to 
the study of the propagation of EM waves in inhomogeneous 
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media. In particular, this method allows one to study geome- 
tries of arbitrary complexity, it can deal with frequency depen- 
dent dielectric functions (metallic inhomogeneous structures) 
in a natural way, discontinuities in the dielectric function are 
not especially detrimental for convergence of the method, and 
the quantities are already calculated in the stationary regime. 
The only shortcoming of the method is its extensive com- 
puter memory usage. However, the demand of this resource is 
quite insensitive to the presence of defects, which could ren- 
der this technique very advantageous for studying disordered 
lattices'^*. Therefore, at the present status of the photonic 
band gap materials field, it is important to assess whether this 
method could be useful to investigate the properties of pho- 
tonic crystals from both the fundamental and applied points 
of view. This is precisely one of the aims of this paper: we 
report a comprehensive study of both the band structure and 
transmission coefficients of various two-dimensional photonic 
crystals using the finite element method. The other main goal 
is to study the feasibility of this method to analyze resonant 
cavities in photonic crystals and their associated quality fac- 
tors. 

The structure of the paper is as follows: in the next sec- 
tion we revisit the equations that describe the propagation 
of electromagnetic fields in two-dimensional inhomogeneous 
media and, in particular, in photonic crystals. Also, the bound- 
ary conditions required for each particular simulation are pre- 



sented. Section III summarizes the results for topologies with 



the symmetry of the square and triangular lattices. Section IV 
examines the dispersion relations for crystals with point de- 
fects for both the square and triangular lattices. Then, defect 
mode tunability is briefly discussed. Section [V] addresses fi- 
nite extent point defect PCs wherein issues such as the effect 
of supercell dimensionality in the formation of photonic local- 
ized states and, hence, the improvement of its quality factor 



are discussed. Finally, in section VII we state the conclusions 
of this work. 



II. THE MATHEMATICS BEHIND THE PROPAGATION 
OF EM WAVES IN INHOMOGENEOUS MEDIA 

The propagation of an electromagnetic wave through an in- 
homogeneous medium (such as the photonic crystals consid- 
ered in this chapter) is described by Maxwell's equations. In 
particular, for two-dimensional media, any electric or mag- 
netic field can always be expressed as a linear combination 
of a transversal electric (TE) and a transversal magnetic fields 
(TM). In the first case, the electric field is perpendicular to 
the photonic crystal plane -whereas the magnetic field is con- 
strained into this plane- and the sourceless Maxwell's equa- 
tions are reduced to a Hemholtz equation for the electric field 
given by 



"^^E-if) + kler(r)E,(r) = 0, 



(1) 



where E~(r) is the z-th component of the electric field at posi- 
tion r, erir) is the inhomogeneous relative dielectric constant 
of the photonic crystal, and = co/c with w the angular fre- 
quency of the incident electric field and c the speed of light 



in free space. In writing down equation ([T|l, it has been as- 
sumed that the photonic crystal is non-magnetic (/v,- = 1) and 
non-conducting (cr = 0). Once equation ([T]) is solved, the 
time-harmonic electric and magnetic fields are easily calcu- 
lated as 



E(f,t)^E,{i)e-"^'z 



H(f,f) = — Vxe,(r)£'(r,f). 



(2) 
(3) 



On the other hand, in the case of TM polarization, the mag- 
netic field is perpendicular to the photonic crystal plane - 
whereas the electric field is constrained into this plane- and 
the sourceless Maxwell's equations reduce to a Hemholz 
equation for the magnetic field given by 

where HAf) is the z-th component of the magnetic field at 
position r. The time-harmonic electric and magnetic fields 
are easily calculated once equation Q is solved and they are 
given by 



Eitt)^- 



koe(r) 



V X H(r, f). 



(5) 
(6) 



An important aspect of any electromagnetic simulation is the 
use of appropriate boundary conditions at the interfaces. On 
the one hand, for the photonic band structure calculations re- 
ported below, it is necessary to implement boundary condi- 
tions that mimic an infinite simulation domain together with 
the periodicity of the photonic crystal lattice. An infinite 
medium is simulated by using perfect magnetic (PMC) or per- 
fect electric conductor (PEC) boundary conditions that mirror 
the simulation domain. The former condition is used for TE 
polarization of the EM field and ensures that the component of 
the magnetic field tangent to the boundary is identically zero 
at the outer interface, that is. 



(7) 



where n is a unit vector perpendicular to the outer simulation 
domain surface at each point. The later is used for TM polar- 
ization of the EM field and ensures that the component of the 
electric field tangent to the boundary is identically zero at the 
outer interface, that is. 



n X 



(8) 



The periodicity of the photonic crystal lattice is ensured by 
an adequate use of Bloch's theorem at the boundaries of the 
photonic crystal unit cell. This theorem states that when the 
electric (magnetic) field propagates from one point in the PC 
to another one separated from the previous one by a lattice 
vector, ^, the only effect on the EM field is a change of its 
phase. 



E,(r + R) = e'^^E.ir) 



(9) 
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and 

H,(r + R) = e''''^H,(r) (10) 

for TE and TM polarization, respectively. In these expres- 
sions, R is a vector of the photonic crystal lattice and k is the 
wavevector of the electromagnetic wave. On the other hand, 
for the transmittance calculations reported below, finite clus- 
ters in the direction parallel to the incident wave vector were 
used, whereas the clusters are of infinite extension in the per- 
pendicular direction. In order to mimic such a material, PMC 
and PEC boundary conditions were used for the outer inter- 
faces that limit the simulation domain in the direction perpen- 
dicular to the incident wave vector for TE and TM polariza- 
tion, respectively. For the interfaces at which the EM wave en- 
ters and exits the cluster, the situation is a little more involved 
because it is necessary to avoid unphysical reflections due to 
the finite size of the cluster and, thus, perfectly matched lay- 
ers were used at these interfaces. The equations that describe 
such boundaries are given by 

z ■ n X (V X £,£) - i/3E: = -2i/3Eo, (1 1) 

z ■ « X (V X H,z) - i/3H, = -UpHa-, (12) 

where isoz and //q; are the initial values of the electric and 
magnetic fields at the boundaries, respectively, and fi - is 
the propagation constant. The upper condition applies to TE 
polarization, whereas the lower one corresponds to TM polar- 
ization. If the electric field is an eigenmode of the boundary, 
the boundary is exactly non-reflecting. 

III. SIMULATION RESULTS FOR PHOTONIC CRYSTALS 
BASED ON SQUARE AND TRIANGULAR LATTICES 

The first structure analyzed in this work is a photonic crys- 
tal made of dielectric circles whose centers occupy the posi- 
tions of a square lattice and is depicted in the inset of Fig. [T| 
The dielectric material was assumed to be linear, isotropic, 
and non-magnetic' . The dielectric constant of the rods has a 
value of 9. The ratio ^, where r is the radius of the cylin- 
ders and a the lattice parameter, was taken as 0.38. This 
well known structure was first investigated by McCall and 
coworkers in order to compare the predictions of theory with 
experimental results with regards to the localization of light 
in strongly scattering media. In the present context, we have 
studied this topology in order to check how well the EE 
method fares when compared with the plane wave method 
that, as stated above, is the one commonly used to calculate 
photonic band structures. The first nine photonic bands were 
calculated for transversal electric (TE) polarization along the 
path that delimits the irreducible part of the 1 st Brillouin Zone 
(IBZ). For the EE calculation, the square unit cell was di- 
vided in 3720 mesh elements and periodic boundary condi- 
tions given by Bloch's theorem were implemented. For the 
MPB calculation, a resolution of 64 x 64 (= 4096) grid ele- 
ments were used and the dielectric constant was average over 
9 grid points. The resulting photonic band structure is de- 
picted in Fig.[T] As it is customary, the dimensionless quantity 



(jjajlnc - a/A has been used to characterize the frequency of 
the incident EM wave, where co is the frequency of the inci- 
dent EM wave and A the associated wavelength. The corre- 
sponding eigenmodes of the z-component of the electric field, 
E~, were also calculated at the F, M, and K points of the IBZ 
and they are shown also in figure [T] It is clear from that figure 
that the band structure calculated with the FE method faith- 
fully reproduces the one calculated with MPB to its minimum 
details. There are three photonic gaps in the band structure of 
this lattice whose sizes coincide with the calculated ones with 
MPB. Also, the modes calculated with the FE method closely 
resemble those calculated with MPB up to a trivial symme- 
try operation or linear combination of degenerate modes. In 
addition, a quantity that can be readily calculated with the 
FEM method is the transmittance of a finite photonic crys- 
tal. The transmittance of the cluster is calculated using the 
usual approach of power integration along a domain bound- 
ary. The boundary conditions for the simulation domain were 
set as follows: on the input and output boundaries (left and 
right boundaries) a perfectly matched layer was used to avoid 
spurious reflections from non-physical boundaries. The z- 
component of the electric field was set to 1 and at the initial 
time of the simulation on the input and output boundaries, re- 
spectively. On the other hand, perfectly magnetic conductor 
boundary conditions that set the tangential magnetic compo- 
nent of the magnetic field to zero were used for the upper and 
lower boundaries" in order to mimic an infinite stripe in that 
direction. 

Another important class of two-dimensional photonic crys- 
tals is the one formed by those structures based on the triangu- 
lar lattice because this is the symmetry most commonly used 
for practical application^. For this reason it is important to 
test whether the FE method produces correct results for this 
non-orthogonal lattice. A photonic crystal made of dielec- 
tric cylinders whose centers occupy the sites of an hexagonal 
lattice of lattice constant a (see the insets in Fig. [2|i was an- 
alyzed. A value of 12 was used for the dielectric constant 
of the non-empty part of the photonic crystal. The radius of 
the cylinders is 0.12a. The simulation setup was very similar 
to the one used for the square lattice: For the FE band cal- 
culation, the unit cell was divided into 1086 mesh elements 
and periodic boundary conditions were implemented across 
the boundaries. For the transmittance calculation in the FM 
direction, a cluster comprising 12 x 8 unit cells (dielectric 
rods) was divided into 60748 mesh elements. For the trans- 
mittance calculation in the TK direction, a cluster comprising 
11x10 unit cells was divided into 67282 mesh elements. The 
same boundary conditions as for the square lattice calcula- 
tions described above were used. The photonic band structure, 
transmittances in the FM and FK directions, and E. field pat- 
terns for TE-polarized EM waves are reported in Fig. [2] The 
band structure calculated with the FE method faithfully repro- 
duces the one obtained with MPB. This structure possesses 
one gap between the first and second bands for TE-polarized 
EM waves that can be clearly seen in the transmission spec- 
tra for light propagating in both the FM and FK directions as 
wide opaque regions. There are also some partial gaps, as the 
ones occurring between the second and third bands and the 
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FIG. 1: (color online). On the left: band structure calculated with the MPB (blue circles) and COMSOL (red line) software 
packages for TE-polarized EM waves. On the right: E, patterns in the unit cell calculated with MPB (upper row in each 
rectangle) and COMSOL at the special symmetry points F, X, and M. 



fourth and fifth bands in the TM direction that, however, are 
not present in the TK direction. Also, there are some opaque 
regions not related to gaps on the band structure but rather to 
the existence of uncoupled modes, such as the one associated 
to the 5th band, that is uncoupled in the VK direction. 

Even though that is not the main purpose of this paper, it is 
illustrative to compare to some extent the accuracy and speed 
of the EE calculations with the ones performed by using the 
PW method. To estimate the relative accuracy, we computed 
the percent error in the eigenvalue calculated at the X point for 
the ninth band for different discretizations of the square unit 
cell. In the FE case, meshes with 254, 414, 928, and 1502 
elements were used, whereas for the calculations done with 
MPB, resolutions of 16 x 16 (= 256) grid elements, 32 x 32 
{- 1024) grid elements, 48 x 48 (= 2304) grid elements, and 
64 X 64 {- 4096) grid elements were used. We took the eigen- 
value calculated with MPB by using a resolution of 256 x 256 
and a mesh size of 25 as the exact one. The result of this com- 
parison can be seen in Fig. [3^. It is noticeable that the FE 
method gets a better accuracy with coarser discretizations of 
the lattice than the PW does. This is due to the use of second 
order Lagrange elements. However, the differences between 
both methods should be negligible for most applications. Of 
course, this better accuracy comes at a price and the simula- 
tion times for the FE calculations are longer by a noticeable 
factor than those for the MPB ones, as can be seen in Fig.[3]3, 
where we have depicted the evolution of the simulation run- 
time with the number of mesh elements and grid elements for 



the FE and MPB calculations, respectively. 



IV. DEFECT STATES IN 2D PHOTONIC CRYSTAL 
SUPERLATTICES 

In complete analogy with semiconductors, the physics of 
disordered photonic crystals is even richer and their potential 
applications even more promising. Among these, photonic 
crystals in which defects are placed in a controlled way form- 
ing point or line defects are especially important for telecom- 
munication applications. Their importance is rooted on the 
fact that these types of defects can lead to the formation of 
localized states inside the photonic gap that allows one to lo- 
calize light around the defects. 

Point defect photonic crystal configurations have been 
widely studied as promising candidates for the enhancement 
of strong coupling between the resonance cavity and quan- 
tum dots light emitters^-* or to reduce the lasing threshold of 
a certain laser emittei^'. Besides, they have been described 
as key elements for several applications in many areas of 
physics and engineering such as enhancing high directivity 
antennas'*, designing low power consumption and highly tun- 
able optical buffering devices'^ or constructing new PC- 
VCSEL lasers'' as well as for biosensing applications^", just 
to cite some examples. 

In the simulations reported in this section, periodicity of 
the lattice has been intentionally broken and hence a defect is 
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Transmittance (FM) Transmittance (FK) 




FIG. 2: (color online) Left panel: Transmittance for TE-polarized waves propagating in the TM direction of the hexagonal 
photonic crystal cluster depicted in the inset. Central panel: Comparison between the band structure calculated with COMSOL 
(red line) and MPB (blue circles) for TE-polarized EM waves along the boundary of the irreducible part of the IBZ. The inset 

shows the unit cell used for the calculations. Right panel: Transmittance for TE-polarized waves propagating in the FK 
direction of the hexagonal photonic crystal cluster depicted in the inset. Inside the green rectangle: Comparison between the E- 
field patterns at the M point calculated with MPB (on the left) and FEM (on the right) for the first eight bands (the band index 
increases from bottom to top). A portion of the photonic crystal that contains 3x3 unit cells is displayed in order to show the 

hexagonal symmetry of the modes. 



introduced into the otherwise perfectly ordered dielectric dis- 
tribution, giving rise to localized states within the band gap re- 
gion. When a single rod is removed from the dielectric lattice, 
light bounces back and forth in the disordered area, trapped 
by the surrounding band gap, whilst no other leakage mech- 
anism is present. The resulting photonic crystal can still be 
classified on the basis of unit cell calculations, as long as the 
fundamental periodic cell hosts sufficient rods around the im- 
perfection to make sure that the defect state does not overlap 
with neighboring replicas of itself. Furthermore, in these sin- 
gle rod perturbed cells, rotational symmetry remains invariant 
and thus comparative studies with previous perfect square and 
triangular lattice based calculations can still be performed at 
the edge of the IBZ. 



A. A square lattice of dielectric rods with a point defect: the 
McCall's experiment 

The experimental evidence of spatial mode localization in 
an ordered dielectric lattice of rods in an air background was 
first given by McCall et al.'^. In the present context, we have 
studied this topology by EE method and compared it to the 
predictions of the PWE method**. Figure [4] depicts the result- 
ing dispersion diagram calculated for the experimental setup 
raised by McCall, where a single rod has been suppressed 
amidst the square lattice of dielectric rods, for which the nor- 
malized rod radius is 0.38a. The rod dielectric constant has 
been set to 9, as in the original McCall's work. A 5 x 5 super- 
cell has been used for the FEM calculations reported in this 
section. It has been discretized into 95903 mesh elements. A 
3096 element grid has been used for the corresponding MPB 
calculations. The real part of the z-component of the electric 
field for the defect mode located at a/ A = 0.4686 was calcu- 
lated at the r point of the IBZ. The main drawback of the su- 
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FIG. 3: (color online) (a) Evolution of the relative error with 
the number of elements in which the lattice is discretized. (b) 
Evolution of the simulation run time with the number of 
elements in which the lattice is discretized. 



percell approach is the band-folding effect: redundant bands 
of the unit cell are folded back times (A^ being the linear 
dimension of the supercell). This fact leads to larger compu- 
tational times since the amount of eigenvalues to be solved 
grows ~ A^. We have therefore solved eigenvalues for 211 
points in /T-space for 100 bands by both methods. Experi- 
mental results obtained by McCall and coworkers, MPB cal- 
culations, and FEM results fully agree as it is shown in Figs.|4] 
andlD 

Two modes make their way across the upper band gap re- 
gion but only the defect mode depicted in Fig.|4^ concentrates 
its energy around the missing rod. 



The density of states (DOS) shows that perturbation of a 
single rod induces a series of sharp Dirac peaks centered at the 
frequencies where the defect state occurs. DOS calculations 
have been also computed via FEM by randomly sampling k- 
points constrained to the UTeducible portion of the IBZ for 
each lattice. In the long wavelength limit, this quantity clearly 
exhibits the linear behavior expected for propagation in an ho- 
mogeneous two-dimensional dielectric medium. 



B. A triangular lattice of dielectric rods with a point defect 

An analogous situation can be achieved by removing a sin- 
gle rod in a photonic crystal based on the triangular lattice. 

Fig. 1^ shows the corresponding results for the photonic 
crystal parameters used by Smith et al.^^ in their investigation 
of the defect mode structures in the square and triangular lat- 
tices. Rod radii and dielectric constant parameters of the trian- 
gular lattice supercell are the same as in the orthogonal case of 
Fig.|4^ and, once again, FEM calculation and results obtained 
by means of PWE method coincide. The simulation setup was 
similar to the previously described one for the square lattice, 
but this time discretization of the supercell has been set to 
71376 mesh elements. Furthermore an hexagonal supercell 
has been used for the FEM calculations, as its shape matches 
the reciprocal lattice of a periodical triangular arrangement. 
In both cases, square and triangular lattice point defect super- 
cells, the electric field z-component is well localized around 
the defect neighborhood and rapidly decays to small ampli- 
tudes as one moves further from the defect. Also, the eigen- 
function around the lattice irregularity clearly shows an in- 
herent symmetry. Indeed, both lattices present a monopole 
pattern with a single nodal plane through each dielectric rodH. 
The symmetry of such point defect modes is analyzed in detail 
iiP. 

The confinement of the electric field amplitude in the cavity 
is directly related to the band gap strength. Stronger mode lo- 
calization has been observed for wider gap geometries. This 
is a highly desirable effect since the gap-mid-gap ratio can 
be easily tuned by adequately preselecting rod radii and di- 
electric contrast. For rods of dielectric contrast equal to 9 in 
square and triangular lattices, one can maximize the gap-mid- 
gap ratio for TE polarization by just adjusting the rod radii to 
0.17fl. Band diagrams and DOS calculations for 0.17fl radius 
rods for square and triangular lattice are reported in Fig. 
and Fig. |5]3. There, localized defect bands are strongly con- 
fined by a wider band gap. Moreover, these defect modes 
can be spatially translated into frequency by simply choos- 
ing the appropriate index contrast between lattice rods and 
point defect rocffl With regards to this mode tunability some 
important properties of photonic crystal cavities can be con- 
veniently enhanced. Indeed, incremental design procedures 
have successfully been demonstrated by means of algorithmic 
techniques rather than using intuitive trial and error design 
methods, where the defect dielectric function is tuned accord- 
ing to design requirement^^SHZl These facts will be discussed 
elsewhere^^. 
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FIG. 4: (color online) a) On the left: Band structure of a 5 x 5 square lattice periodic supercell calculated with MPB (blue 
circles) and FEM (red lines), wherein a defect state has been excited around ^ - 0.466 ± 0.002. On the right: density of 
photonic states (DOS), calculated by FEM. The band gap region is clearly distinguishable and a weakly localized defect mode 
merging from the upper band can be seen. In the inset, E- patterns for a 5 x 5 supercell calculated by both methods are shown 
for the M symmetry point, b) On the left: band structure of a 5 x 5 triangular lattice calculated with MPB (blue circles) and 
FEM (red line). A single defect state merges st ^ - 0.285 ± 0.003. On the right: density of photonic states, calculated by 
FEM. In the inset, patterns for the F symmetry point. An hexagonal supercell has been used for the FEM calculation . 



V. DEFECT STATES IN FINITE 2D PHOTONIC CRYSTAL 
CLUSTERS 

All the structures studied so far had an infinite extension 
in all space directions. Nevertheless, an engineered photonic 
crystal must be finite sized. This fact limits defect mode con- 
finement factor as stored energy will be finally radiated out- 
side the cavity. This circumstance gives rise to a quantity of 
very much practical importance: the quality factor of a pho- 
tonic crystal resonator, Q, which is just the ratio between the 
stored energy in the photonic crystal and the radiated energy 
per cycle. In the present context it is unnecessary to consider 
the radiated energy in the off-axis direction because the crystal 
has still infinite length in this direction, but inasmuch as the 
in-plane propagation will be limited to periods of dielec- 
tric rods, electromagnetic fields will still be leaked outwards. 
Quantitatively, these energy decay mechanisms are uncorre- 
cted to each other, and so they can be studied separately. On 
the other hand, in the present calculations no losses, absorp- 
tions, nor imperfections have been taken into account and, 
consequently, the quality factor is only limited by the inher- 
ent energy leakage from the cavity to the radiative modes. 

Transmittance diagrams have been obtained using the FEM 



but, in contrast with the transmittances reported above for 
ideal lattices, in this section a finite length design is regarded. 
Therefore, PMC and PEC boundary conditions were substi- 
tuted by PML interfaces in the direction perpendicular to the 
incident wave vector as well as at the interfaces at which the 
EM wave enters and exits the structure. Accordingly, trans- 
parent influx conditions were imposed to the outer boundaries 
of the finite cluster and a TE polarized source propagating in 
the r - X direction was considered. Spurious reflections are 
thus prevented from being injected back into the simulation 
domain by locating these artificial perfectly matched layers 
all around the simulation domain boundaries. The equations 
that describe such boundaries are given by 

z ■ n X (V X E,z) - i koE. = -ika{\ - h ■ «)£oze"*"° '' (13) 

z ■ n X (V X H,z) - ikoH, = -ikoil - ko ■ n)i/o;e"*° '', (14) 

where Eq, and Hq, are the initial values of the electric and 
magnetic fields at the boundaries, respectively, and is the 
propagation constant. A set of monochromatic plane waves 
is excited for each frequency in the simulation domain. This 
way, a point defect is introduced in the topologies character- 
ized in previous sections. The resulting transmittance results 
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FIG. 5: (color online) a) On the left: band structure of a 5 x 5 square lattice periodic supercell calculated with MPB (blue 
circles) and FEM (red lines). The rod radii has been set to O.lTa, which ensures a better localization of the unique defect mode. 

On the right: density of photonic states, calculated by FEM. The defect mode clearly shows up in the gap region, since the 
defect state is strongly confined. In the inset, the patterns for a 5 x 5 supercell calculated with both methods are shown for 
the M symmetry point. All these patterns match the ones presented in Fig. [4^ but the modal volume decreases significantly for 
0.17fl rod lattices, b) Triangular lattice made of Q.lla rods where the central rod has been removed. The band gap increases and 

so it does the defect mode localization around the point defect. 



are depicted in Fig.|6] 

When the size of the cluster is enlarged, by adding subse- 
quent periods to the lattice, the transmittance spectra clearly 
reveals a sharp peak centered at the defect mode frequency. 
When one considers a larger cluster, the defect mode gets 
narrower, i.e. it converges to the previously discussed peri- 
odic case (Fig.|6]). Additionally, the defect mode frequency is 
shifted due to the overall dielectric percentage change in the 
finite cluster which is indeed very convenient since it allows 
one to tune the resonance frequency according to dielectric 
index changes. 

For the particular case of a square arrangement of dielectric 
rods with radius 0.38fl and dielectric constant of 9, wherein its 
dimensions have been gradually modified, the energy stored 
in the defect mode centered around | - 0.4686 is exchanged 
with the upper band energy situated around ^ = 0.49. The 
quality factor has been calculated by means of the transmit- 
tance response for each crystal, since by definition 

fc 



FWHM 



(15) 



where represents the central frequency of the defect mode 
and FWHM the difference between the two values in the trans- 



mittance function at which the transmittance is equal to the 
half of defect mode maximum amplitude. For each quality 
factor calculation, the FWHM has been accurately determined 
by means of the Brent's algorithm^', which combines bisec- 
tion, regula falsi, and inverse quadratic interpolation methods 
for root finding. In contrast with the accepted wisdom, this 
method has proved to be highly reliable and more efficient 
than using FDTD, as we explain below in more detail. As 
the amount of dielectric rods surrounding the defect increases, 
the rate of the energy loss within the cavity relative to the en- 
ergy confined in it decays exponentially. Figure [7] depicts the 
quality factor increment for these square and triangular clus- 
ters. According to that figure, the defect mode bandwidth de- 
creases exponentially with the addition of new periods due to 
the strengthening of the band gap effect, which is the main 
contribution to the enhancement of the quality factor. The on- 
set of the leakage mechanism transforms the localized modes 
transmittance spectra into distorted Lorentzian peaks. These 
peaks progressively tend to a Lorentzian curve when the de- 
fect mode gets closer to the mid-gap frequency and its FWHM 
gets narrower whenever increases. 



9 




0.44 0.46 

wa/2nc 



0.38 0.42 
(Dfl/27IC 




aa/2nc (aa/2nc 

FIG. 6: (color online) a) Top: transmittance diagrams for different cluster dimensions in a rectangular arrangement of dielectric 
rods with a central defect. On the left, rod radii is set to 0.38fl while on the right ^ is Q.lla. In the inset, detailed defect mode 
curves are depicted.(b) Down: analogous transmittance diagrams in a triangular cluster 



VI. TIME DOMAIN APPROACH: FDTD 

As FDTD is the most popular method for the theoretical 
description of light propagation in these systems, we used 
FDTD in order to assess and compare the Q factor calculations 
attained so far using FEM. For the present manuscript, the 
FDTD simulations were performed using a freely 

available implementation package of this method. The suc- 
cess of FDTD method is due to its flexibility and to its adapt- 
ability to irregular or aperiodic geommetries. More generally, 
the FDTD method can compute a large number of frequencies 
at once and even extract modes of the spectrum. 

In FDTD space and time is divided into a finite rectangular 
grid and then fields are evolved in time using discrete steps. 
However, FDTD has serious limitations when the computa- 
tional domain is finite whilst FEM convergence time is rather 
insensitive to this fact. In FDTD the computational domain 
must be terminated with some boundary conditions as it is the 
case in the FEM and, in order to simulate open boundaries 
in finite clusters, a wave absorbing mechanism, such as the 
split-field PML proposed by Berenge^has been used. Such 
an artificial region is needed so as to absorb outgoing spuri- 
ous waves from the computational grid rather than reflecting 
them back into the photonic crystal. Moreover, quality fac- 
tor calculations are very sensitive to the size of the computa- 



tional grid and, thus, if the injected pulse experiences spuri- 
ous reflections from the domain boundaries, radiated normal- 
incident waves will not be the dominant ones and they will 
be mixed with reflected waves. MEEP implements Uniax- 
ial Perfectly Matched Layers (UPML)^^ in order to absorb 
outgoing spurious waves from the computational. Along this 

artificial medium, the PML is expressed in terms of effective 
1 Its' 

anisotropic e and v . 

For the present FDTD computation, a point Gaussian 
source has been placed inside the cavity. This excitation must 
be short enough (broad bandwidth) to excite the defect mode 
for each cluster When this Gaussian source is switched on, 
the field grows and after some time this source is extinguished. 
Subsequently, a resonance effect occurs and the electromag- 
netic fields bounce back and forth for a limited amount of 
optical periods. Meanwhile, the energy trapped around the 
defect exhibits an exponential time decay (see Fig.[8]l. At this 
point, if one takes into account the slowly varying component 
only, i.e., the envelope of the electric field norm, the quality 
factor is determined by the ratio of the stored power divided 
by the loss power after one cycle: 
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FIG. 7; (color online) Top: quality factor estimation for a number of different sized square clusters of 0.38a radii dielectric rods 
(left) and Q.lla rods (right). Down: analogous computation results obtained for triangular arrangements. 



In this type of calculations, the proximity of the UPML to 
the cluster significantly determines the decaying factor of the 
existing resonance. In the FDTD calculations reported in Ta- 
bleUan UPML of thickness 2a surrounds each structure. Be- 
sides, the distance between the UPML layer and the edge of 
the computational grid must be adequately tuned in order to 
avoid an unphysical field decay. According to this, a distance 
of 0.8 < d < 1.1 has been used. Sampling the electric field 
response with a period relative to the resonance frequency, 
results in a quasi-periodic step function that induces some un- 



certainties in the Q factor determination if ( 16 1 is directly ap- 
plied. These fluctuations could be due to the abruptly broken 
translational symmetry of the clusters. However, the marked 
decaying behavior of the electric field data can be easily fil- 



tered and interpolated. Therefore, by iteratively applying ( 16 1 



to the interpolated sample field, an average Q factor and the 
corresponding theoretical error is obtained. Finally, in order 
to compare with previous Q factor calculations obtained via 
FEM transmittance results and FDTD decaying field value re- 
lations, an harmonic inversion of time signals (Harminv)^"^ 
implementation has been used. This software package uses 
a filter diagonalization method (FDM) to find the deconvolu- 
tion of sinusoidal functions near a given frequency interval. 
The simulation setup is the same as the one used in the stan- 
dard FDTD calculations: a broad Gaussian point source lo- 
cated in the defect excites TE modes in the cavity and PML 
layers surrounding the cluster. In particular, a | = 0.02 pulse 
width Gaussian source centered at the resonant frequency of 
each structure has been excited. When the source vanishes, 
Harminv performs the signal processing of the fields in the 



cavity. This way, it identifies the frequencies and decay rates 
of the excited resonant modes. This method is a fast and accu- 
rate way to determine the Q factors but still requires waiting 
until the fields have evolved and decayed to a certain small 
value. Computation time scales also with the cluster size, but 
this fact can often be overcome if symmetries are used to re- 
duce the size of the computational cell. Besides, choosing a 
broad Gaussian source can yield to the appearance of spuri- 
ous solutions but, at the same time, the source must be broad 
enough to excite the resonant frequency. 

It is important to stress the remarkable agreement between 
FEM results, FDTD and Harminv calculations of the Q fac- 
tor Table [l] summarizes the results obtained using the three 
methods described in this report. Two dimensional FEM cal- 
culation of Q factors offers a very satisfactory agreement with 
FDTD and Harminv results whilst providing substantial gain 
in solution robustness and efficiency. Indeed, FDTD is a pow- 
erful tool to calculate resonant frequencies and quality fac- 
tors for complex cavity structures. However, it is very ineffi- 
cient because one must discard many simulation cycles before 
reaching the stationary regime. Also, there are many heuris- 
tic factors that enter into the FDTD simulation process, such 
as the thickness of the UPML, the excitation source location 
and both its central frequency and width, that make the results 
from this method not as reliable as one would want. Moreover, 
when computing resonant modes in time domain especial care 
must be taken to avoid choosing an excitation source nearly 
orthogonal to the resonant mode because FDTD is likely to 
miss it, otherwise. In this regard, FEM can be seen as an effi- 
cient, reliable and more rigorous alternative to FDTD for the 
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FIG. 8: (color online) Top: evolution of the electromagnetic field inside the defect rod for a 0.38a arrangement and for 
different cluster dimensions. A point source is excited in the cavity but after some periods the source is extinguished and 
although the electromagnetic field still remains in the cavity it experiences an exponential decay. Down: analogous 

computation results obtained for triangular defect cluster 



analysis of quality factors and resonant modes in complex di- 
electric material. 



VII. CONCLUSIONS 

The present manuscript reports a comprehensive study of 
the photonic properties of several two-dimensional photonic 
crystals and finite clusters by using the finite element method. 
The main result coming out from these calculations is that 
the FEM allows one to reproduce the results obtained with 
the well-known plane wave and FDTD methods, but it has 
many advantages not present in the others. In contrast with 
frequency-space based approaches, FEM also deals in a nat- 
ural way with finite clusters and aperiodic structures of arbi- 
trary complexity. To prove this, we have calculated the band 
structure of periodic photonic crystals based on the square and 
triangular lattices. It is found that the band structures calcu- 
lated in this way are almost indistinguishable of those calcu- 
lated with the well known MPB package. Also, the modes 
calculated with FEM closely resemble those calculated with 
the PWE method. It is noticeable that using a coarser dis- 
cretization FEM results are more acurate than the ones given 



by PWE. 

Moreover, the transmission coefficients of a number of fi- 
nite clusters of the aforementioned lattices were calculated 
along different directions in reciprocal space and for both TE 
and TM polarizations. The features (such as the position and 
width of the photonic gaps) of these transmittances agree quite 
well with the band structures and once again the accuracy of 
the band and mode calculations is very good when compared 
with the PWE method. Therefore, these results demonstrate 
that the FEM method can be a very useful general purpose 
method for investigating photonic crystals. This fact is further 
stressed by the results obtained for PCs containing a single 
defect, where DOS and dispersion diagrams have been cal- 
culated. There, experimental results obtained by McCall and 
coworkers, MPB calculations, and FEM results fully agree. 
As seen in these cavities, the confinement of the electric field 
amplitude strongly depends on the band gap strength of the 
underlying structure and, thus, wider gap geometries support 
higher mode localization. In addition, former calculations 
have been reproduced for finite length point defect clusters. In 
this context, the localization of light around the defect region 
has been quantified by accurately determining the quality fac- 
tor using FEM, FDTD, and Harminv procedures for different 
cluster dimensions. 
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Square lattice r=0.38a Square lattice r=0. 17a 
Cluster size FEM FDTD Harminv FEM FDTD Harminv 
5 72 80+1 88 108 113±13 107 

7 376 350±3 347 722 727±6 723 

9 1800 1773±24 1765 4310 4312±31 4308 
11 8350 7677±83 7674 24500 24561 ±74 24559 

Triangular lattice r=0.38fl Triangular lattice r=0.17a 
5 97 107±2 103 134 130+1 126 

7 358 320±1 316 1280 1171±8 1166 

9 1240 121 1±5 1209 11100 11241±42 11238 
11 5270 5656±74 5650 90200 83903±84 83892 



TABLE I: Quality factor calculations for (liferent cluster sizes in square and triangular lattice rod-type PC single-defect cavities 

using FEM, FDTD, and harmonic inversion methods. 



FEM is proven to be an effective and stable tool for point 
defect cluster quality factor calculation, wherein for each 
quality factor calculation, the FWHM has been accurately de- 
termined by means of Brent's algorithm. This technique per- 
mits to determine the essential information needed for Q fac- 
tor calculation in a speedy and computationally effective way. 
In addition, the leakage mechanism of PC cavities transforms 
the transmittance spectra into an almost Lorentzian peak and, 
therefore, with few transmittance calculations one can obtain 
the entire point defect transmittance response. It is notewor- 
thy, with regards to the the point defect PCs addressed above, 
the fact that the three methods accurately reproduce the Q fac- 
tor for each topology. However, FDTD based methods have 
serious drawbacks. On the one hand, the width and the prox- 
imity of the absorbing layers significantly determine the de- 
caying behavior of the trapped fields inside the cavity, and so, 
the parameters of the PMLs must be carefully chosen for each 
structure. In the aforementioned point defect structures, the 
source must be placed inside the cavity, which in some cases 
can be an unrealistic situation. Furthermore, the bandwidth of 
the source must be set intuitively so as to excite only the defect 
mode. On the other hand, we believe that among these meth- 



ods FEM is desirable for the calculation of Q factor due to its 
high numerical efficiency and stability because in FDTD, after 
the source is extinguished, one must wait an uncertain amount 
of time until the fields evolve and decay. In fact, FDTD gives 
quite accurate values for both the resonant frequency and the 
Q factor but, for higher Q values, the slope of the electromag- 
netic field inside the defect is nearly zero and hence, conver- 
gence time and and numerical errors increase drastically. 

To summarize, the results presented in this manuscript 
demonstrate that the finite element method is an effective, sta- 
ble, robust, and rigorous tool to study light propagation and 
confinement in both periodic and aperiodic dielectric photonic 
crystals. Furthermore, we expect that these advantages can be 
extrapolated to systems in which the optical constants are fre- 
quency dependent, such as hybrid metallo-dielectric photonic 
crystals. 
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